WGCNA Top Modules Exploratory Analysis

library(tidyverse) ; library(reshape2) ; library(glue) ; library(plotly) ; library(dendextend)
library(RColorBrewer) ; library(viridis) ; require(gridExtra) ; library(GGally) ; library(ggpubr)
library(expss)
library(polycor)
library(foreach) ; library(doParallel)
library(knitr) ; library(kableExtra)
library(biomaRt)
library(clusterProfiler) ; library(ReactomePA) ; library(DOSE) ; library(org.Hs.eg.db)
library(WGCNA)

SFARI_colour_hue = function(r) {
  pal = c('#FF7631','#FFB100','#E8E328','#8CC83F','#62CCA6','#59B9C9','#b3b3b3','#808080','gray','#d9d9d9')[r]
}

Load preprocessed dataset (preprocessing code in 01_data_preprocessing.Rmd) and clustering (pipeline in 05_WGCNA.Rmd)

# Gandal dataset
load('./../Data/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
rownames(datExpr) = datGenes$ensembl_gene_id
datMeta = datMeta %>% mutate(ID = paste0('X',description))


# GO Neuronal annotations: regex 'neuron' in GO functional annotations and label the genes that make a match as neuronal
GO_annotations = read.csv('./../Data/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>% 
              mutate('ID'=as.character(ensembl_gene_id)) %>% 
              dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
              mutate('Neuronal'=1)


# SFARI Genes
SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_01-03-2020_w_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]

# Old SFARI Genes
old_SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_08-29-2019_w_ensembl_IDs.csv')
old_SFARI_genes = old_SFARI_genes[!duplicated(old_SFARI_genes$ID) & !is.na(old_SFARI_genes$ID),]


# Clusterings
clusterings = read_csv('./../Data/clusters.csv')


# Update DE_info with SFARI and Neuronal information
genes_info = DE_info %>% mutate('ID'=datGenes$ensembl_gene_id, padj = adj.P.Val, log2FoldChange = logFC) %>% 
             left_join(SFARI_genes, by='ID') %>% 
             mutate(`gene-score`=ifelse(is.na(`gene-score`), 'Others', `gene-score`)) %>%
             left_join(old_SFARI_genes %>% rename(`gene-score` = 'old-gene-score') %>% 
                       dplyr::select(ID, `old-gene-score`), by = 'ID') %>%
             mutate(`old-gene-score`=ifelse(is.na(`old-gene-score`), 'Others', `old-gene-score`)) %>%
             left_join(GO_neuronal, by='ID') %>% left_join(clusterings, by='ID') %>%
             mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
             mutate(gene.score=ifelse(`gene-score`=='Others' & Neuronal==1, 'Neuronal', `gene-score`), 
                    old.gene.score=ifelse(`old-gene-score`=='Others' & Neuronal==1,'Neuronal',`old-gene-score`), 
                    significant=padj<0.05 & !is.na(padj))

# Add gene symbol
getinfo = c('ensembl_gene_id','external_gene_id')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl',
               host='feb2014.archive.ensembl.org') ## Gencode v19
gene_names = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), values=genes_info$ID, mart=mart)

genes_info = genes_info %>% left_join(gene_names, by=c('ID'='ensembl_gene_id'))


clustering_selected = 'DynamicHybrid'
genes_info$Module = genes_info[,clustering_selected]

dataset = read.csv(paste0('./../Data/dataset_', clustering_selected, '.csv'))
dataset$Module = dataset[,clustering_selected]
dataset$gene.score = as.character(dataset$gene.score)
dataset$gene.score[dataset$gene.score=='None'] = 'Others' 


rm(DE_info, GO_annotations, clusterings, getinfo, mart, efit)

Selecting Top Modules

plot_data = dataset %>% dplyr::select(Module, MTcor) %>% distinct %>% 
            mutate(alpha = ifelse(abs(MTcor)>0.85, 1, 0.3))

top_modules = plot_data %>% arrange(desc(MTcor)) %>% filter(abs(MTcor)>0.85) %>% pull(Module) %>% as.character

Selecting Modules with a Module-Diagnosis correlation magnitude larger than 0.85 (instead of 0.9 as we did with Gandal’s dataset because the relation is not as strong in this dataset)

The 2 modules that fulfill this condition are #79AF00, #00B813

ggplotly(plot_data %>% ggplot(aes(reorder(Module, -MTcor), MTcor)) + 
         geom_bar(stat='identity', fill = plot_data$Module, alpha = plot_data$alpha) + 
         geom_hline(yintercept =c(0.85, -0.85), color = 'gray', linetype = 'dotted') + 
         xlab('Modules')+ ylab('Module-Diagnosis Correlation') + theme_minimal() + 
         theme(axis.text.x = element_text(angle = 90, hjust = 1)))

The modules consist mainly of points with high (absolute) values in PC2 (which we know is related to LFC), so this result is consistent with the high correlation between Module and Diagnosis, although some of the points with the highest PC2 values do not belong to these top modules

The genes belonging to the modules with the positive Module-Diagnosis correlation have positive LFC values and the genes belonging to the modules with the negative Module-Diagnosis correlation have negative values

pca = datExpr %>% prcomp

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% 
            left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>%
            dplyr::select(ID, external_gene_id, PC1, PC2, Module, gene.score) %>%
            mutate(ImportantModules = ifelse(Module %in% top_modules, as.character(Module), 'Others')) %>%
            mutate(color = ifelse(ImportantModules=='Others','gray',ImportantModules),
                   alpha = ifelse(ImportantModules=='Others', 0.1, 0.6),
                   gene_id = paste0(ID, ' (', external_gene_id, ')')) %>%
            apply_labels(ImportantModules = 'Top Modules')

cro(plot_data$ImportantModules)

	#Total
Top Modules
#00B813	708
#79AF00	399
Others	14351
#Total cases	15458

ggplotly(plot_data %>% ggplot(aes(PC1, PC2, color=ImportantModules)) + 
         geom_point(alpha=plot_data$alpha, color=plot_data$color, aes(ID=gene_id)) + theme_minimal() +
         xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],2),'%)')) +
         ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],2),'%)')) +
         ggtitle('Genes belonging to the Modules with the strongest relation to ASD'))

rm(pca)

---

SFARI Genes

List of top 20 SFARI Genes in top modules ordered by SFARI score and Gene Significance

list_top_SFARI_genes = function(table_data, module) {
  
  t = table_data %>% filter(Module == module & `SFARI score` %in% c(1,2,3)) %>% 
      slice_head(n=20) %>% dplyr::select(-Module, -`Ensembl ID`)
 
  return(t)
}


table_data = dataset %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>%
             dplyr::select(ID, external_gene_id, GS, gene.score, Module) %>% 
             arrange(gene.score, desc(abs(GS))) %>%
             dplyr::rename('Ensembl ID' = ID, 'Gene Symbol' = external_gene_id, 
                    'SFARI score' = gene.score, 'Gene Significance' = GS)

kable(list_top_SFARI_genes(table_data, top_modules[1]), 
      caption=paste0('Top SFARI Genes for Module ', top_modules[1])) %>%
      kable_styling(full_width = F)

Top SFARI Genes for Module #79AF00

Gene Symbol	Gene Significance	SFARI score
CASK	0.7991471	1
DLG4	0.6492071	1
SRCAP	0.3073226	1
MECP2	0.0496312	1
DEAF1	0.2414273	2
TAF6	0.2175760	2
GABRG3	0.1697446	2
OXTR	0.1531516	2
CDH13	0.1529404	2
SNX5	0.0422536	2
CMIP	0.4673635	3
MUC4	0.4548544	3
MBD3	0.4369314	3
AVPR1B	0.4227676	3
GRM5	0.3342928	3
PSD3	0.3328927	3
HTR3A	0.2924446	3
LRRC1	0.2344388	3
ZC3H11A	0.2164467	3
TTN	0.1221977	3

kable(list_top_SFARI_genes(table_data, top_modules[2]), 
      caption=paste0('Top SFARI Genes for Module ', top_modules[2])) %>%
      kable_styling(full_width = F)

Top SFARI Genes for Module #00B813

Gene Symbol	Gene Significance	SFARI score
SON	-0.3279955	1
SHANK2	-0.2358767	1
ARX	-0.2314617	1
VPS13B	-0.2273199	1
PTPN11	-0.1303289	1
KDM3B	-0.1218145	1
ICA1	-0.6947087	2
GGNBP2	-0.6446193	2
CACNA2D3	-0.6297103	2
TERF2	-0.5561551	2
MTOR	-0.4602444	2
RBFOX1	-0.4470894	2
KCNS3	-0.4198392	2
GALNT8	-0.2812489	2
APBB1	-0.2747579	2
CHRNA7	-0.2721684	2
VIL1	-0.2406558	2
SLC7A3	-0.2062360	2
OTUD7A	-0.1204609	2
APBA2	-0.6365906	3

rm(table_data, list_top_SFARI_genes)

---

Module Eigengenes

Since these modules have the strongest relation to autism, this pattern should be reflected in their model eigengenes, having two different behaviours for the samples corresponding to autism and the ones corresponding to control

In all cases, the Eigengenes separate the behaviour between autism and control samples very clearly (p-value < \(10^{-4}\)). Modules with positive Module-Diagnosis correlation have higher eigengenes in the ASD samples and Modules with a negative correlation, in the Control samples

plot_EGs = function(module){

  plot_data = data.frame('ID' = rownames(MEs), 'MEs' = MEs[,paste0('ME', module)], 
                         'Diagnosis' = datMeta$Diagnosis)

  p = plot_data %>% ggplot(aes(Diagnosis, MEs, fill=Diagnosis)) + 
      geom_boxplot(outlier.colour='#cccccc', outlier.shape='o', outlier.size=3) +
      ggtitle(paste0('Module ', module, '  (MTcor=',round(dataset$MTcor[dataset$Module == module][1],2),')')) + 
      stat_compare_means(method = 't.test', method.args = list(var.equal = FALSE), label = 'p.signif',
                        ref.group = 'ASD') +
      ylab('Module Eigengenes') + theme_minimal() + theme(legend.position='none')
  
  return(p)
}


# Calculate Module Eigengenes
ME_object = datExpr %>% t %>% moduleEigengenes(colors = genes_info$Module)
MEs = orderMEs(ME_object$eigengenes)

p1 = plot_EGs(top_modules[1])
p2 = plot_EGs(top_modules[2])

grid.arrange(p1, p2, nrow=1)

rm(plot_EGs, ME_object, MEs, p1, p2)

---

Identifying representative genes for each Module

In the WGCNA pipeline, the most representative genes of each module are selected based on having a high module membership as well as a high (absolute) gene significance, so I’m going to do the same

SFARI Genes don’t seem to be more representative than the rest of the genes

create_plot = function(module){
  
  plot_data = dataset %>% dplyr::select(ID, paste0('MM.',gsub('#','',module)), GS, gene.score) %>% 
              filter(dataset$Module==module)
  colnames(plot_data)[2] = 'Module'
  
  SFARI_colors = as.numeric(names(table(as.character(plot_data$gene.score)[plot_data$gene.score!='Others'])))
  
  p = ggplotly(plot_data %>% mutate(gene.score = ifelse(gene.score =='Others', 'Not in SFARI', gene.score)) %>% 
               ggplot(aes(Module, GS, color=gene.score)) +
               geom_point(alpha=0.5, aes(ID=ID)) +  xlab('Module Membership') + ylab('Gene Significance') +
               ggtitle(paste0('Module ', module, '  (MTcor = ', 
                              round(dataset$MTcor[dataset$Module == module][1],2),')')) +
               scale_color_manual(values=SFARI_colour_hue(r=c(SFARI_colors,7))) + theme_minimal())
  
  return(p)
}

create_plot(top_modules[1])

create_plot(top_modules[2])

rm(create_plot)

Top 20 genes for each module

Ordered by \(\frac{MM+|GS|}{2}\)

There aren’t that many SFARI genes in the top genes of the modules

create_table = function(module){
  top_genes = dataset %>% left_join(genes_info %>% dplyr::select(ID, external_gene_id), by='ID') %>% 
              dplyr::select(ID, external_gene_id, paste0('MM.',gsub('#','',module)), GS, gene.score) %>% 
              filter(dataset$Module==module) %>% dplyr::rename('MM' = paste0('MM.',gsub('#','',module))) %>% 
              mutate(Relevance = (MM+abs(GS))/2, 
                     gene.score = ifelse(gene.score =='Others', 'Not in SFARI', gene.score)) %>% 
              arrange(by=-Relevance) %>% top_n(20) %>% 
              dplyr::rename('Gene Symbol' = external_gene_id, 'SFARI Score' = gene.score)
  return(top_genes)
}

top_genes = list()
for(i in 1:length(top_modules)) top_genes[[i]] = create_table(top_modules[i])

kable(top_genes[[1]] %>% dplyr::select(-ID), caption=paste0('Top 20 genes for Module ', top_modules[1], 
      '  (MTcor = ', round(dataset$MTcor[dataset$Module == top_modules[1]][1],2),')')) %>% 
      kable_styling(full_width = F)

Top 20 genes for Module #79AF00 (MTcor = 0.87)

Gene Symbol	MM	GS	SFARI Score	Relevance
WWC1	0.7770464	0.8108212	Not in SFARI	0.7939338
CASK	0.7886111	0.7991471	1	0.7938791
CCDC50	0.8122779	0.7306045	Not in SFARI	0.7714412
ABLIM3	0.7420932	0.7196454	Not in SFARI	0.7308693
RIMS4	0.7437251	0.6842736	Not in SFARI	0.7139994
DLG4	0.7218259	0.6492071	1	0.6855165
RAB11FIP4	0.6854802	0.6449110	Not in SFARI	0.6651956
JUN	0.6378685	0.6890684	Not in SFARI	0.6634684
PRR5	0.6026740	0.7175114	Not in SFARI	0.6600927
MXD4	0.6532816	0.6281044	Not in SFARI	0.6406930
RBM39	0.7523953	0.5196978	Not in SFARI	0.6360465
CCDC102A	0.6342699	0.6362624	Not in SFARI	0.6352661
SERTAD1	0.6300838	0.6313406	Not in SFARI	0.6307122
TMEM132A	0.6166367	0.6436432	Not in SFARI	0.6301399
TNFRSF1A	0.7071841	0.5462540	Not in SFARI	0.6267190
SBF2	0.6532819	0.5634664	Not in SFARI	0.6083741
KIAA1522	0.5645191	0.6477544	Not in SFARI	0.6061367
COMMD2	0.5308295	0.6683542	Not in SFARI	0.5995919
KCNK12	0.6686929	0.5259704	Not in SFARI	0.5973317
MSN	0.6179045	0.5757439	Not in SFARI	0.5968242

kable(top_genes[[2]] %>% dplyr::select(-ID), caption=paste0('Top 20 genes for Module ', top_modules[2], 
      '  (MTcor = ', round(dataset$MTcor[dataset$Module == top_modules[2]][1],2),')')) %>% 
      kable_styling(full_width = F)

Top 20 genes for Module #00B813 (MTcor = -0.91)

Gene Symbol	MM	GS	SFARI Score	Relevance
CIRBP	0.8238270	-0.8158820	Not in SFARI	0.8198545
EIF3K	0.7973822	-0.7788495	Not in SFARI	0.7881159
FIBP	0.7705627	-0.8052324	Not in SFARI	0.7878976
UQCRC1	0.8000285	-0.7699530	Not in SFARI	0.7849907
COPG1	0.7859792	-0.7813419	Not in SFARI	0.7836606
PFKP	0.7868055	-0.7437745	Not in SFARI	0.7652900
RTN4IP1	0.8020323	-0.7283473	Not in SFARI	0.7651898
EIF2AK1	0.7984542	-0.7119941	Not in SFARI	0.7552241
PYGB	0.7675212	-0.7128820	Not in SFARI	0.7402016
C17orf49	0.7426747	-0.7237342	Not in SFARI	0.7332044
ICA1	0.7663316	-0.6947087	2	0.7305202
TCEA2	0.7690770	-0.6881259	Not in SFARI	0.7286015
MRPL2	0.7539501	-0.7031087	Not in SFARI	0.7285294
ACTR1B	0.6882590	-0.7683051	Not in SFARI	0.7282821
ALOX12B	0.6670812	-0.7739453	Not in SFARI	0.7205133
OGDHL	0.7048090	-0.7346210	Not in SFARI	0.7197150
PRC1	0.7844209	-0.6470912	Not in SFARI	0.7157560
GLS2	0.7180064	-0.7038931	Not in SFARI	0.7109497
FH	0.6930093	-0.7237622	Not in SFARI	0.7083858
IMMT	0.6590283	-0.7498505	Not in SFARI	0.7044394

rm(create_table, i)

pca = datExpr %>% prcomp

ids = c()
for(tg in top_genes) ids = c(ids, tg$ID)

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% dplyr::select(ID, PC1, PC2, Module, gene.score) %>%
            mutate(color = ifelse(Module %in% top_modules, as.character(Module), 'gray')) %>%
            mutate(alpha = ifelse(color %in% top_modules & ID %in% ids, 1, 0.1))

plot_data %>% ggplot(aes(PC1, PC2)) + geom_point(alpha=plot_data$alpha, color=plot_data$color) + 
              xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],2),'%)')) +
              ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],2),'%)')) +
              theme_minimal() + ggtitle('Most relevant genes for top Modules')

rm(ids, pca, tg, plot_data)

Level of expression by Diagnosis for top genes, ordered by relevance (defined above): There is a visible difference in level of expression between diagnosis groups in all of these genes

create_plot = function(i){
  
  plot_data = datExpr[rownames(datExpr) %in% top_genes[[i]]$ID,] %>% mutate('ID' = rownames(.)) %>% 
              melt(id.vars='ID') %>% mutate(variable = gsub('X','',variable)) %>%
              left_join(top_genes[[i]], by='ID') %>%
              left_join(datMeta %>% dplyr::select(description, Diagnosis),
                        by = c('variable'='description')) %>% arrange(desc(Relevance))
  
  p = ggplotly(plot_data %>% mutate(external_gene_id=factor(`Gene Symbol`, 
                                    levels=unique(plot_data$`Gene Symbol`), ordered=T)) %>%
               ggplot(aes(`Gene Symbol`, value, fill=Diagnosis)) + geom_boxplot() + theme_minimal() +
                      ggtitle(paste0('Top Genes for Module ', top_modules[i], ' (MTcor = ',
                      round(dataset$MTcor[dataset$Module == top_modules[i]][1],2), ')')) + 
                      xlab('') + ylab('Level of Expression') +
                      theme(axis.text.x = element_text(angle = 90, hjust = 1)))
  return(p)
}

create_plot(1)

create_plot(2)

rm(create_plot)

---

Enrichment Analysis

Using the package clusterProfiler. Performing Gene Set Enrichment Analysis (GSEA) and Over Representation Analysis (ORA) using the following datasets:

• Gene Ontology

• Disease Ontology

• Disease Gene Network

• KEGG

• REACTOME

file_name = './../Data/GSEA.RData'

if(file.exists(file_name)){
  load(file_name)
} else {
  ##############################################################################
  # PREPARE DATASET
  
  # Create dataset with top modules membership and removing the genes without an assigned module
  EA_dataset = data.frame('ensembl_gene_id' = genes_info$ID, module = genes_info$Module)  %>% 
               filter(genes_info$Module!='gray')
  
  # Assign Entrez Gene Id to each gene
  getinfo = c('ensembl_gene_id','entrezgene')
  mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl', 
                 host='feb2014.archive.ensembl.org')
  biomart_output = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), 
                         values=genes_info$ID[genes_info$Module!='gray'], mart=mart)
  
  EA_dataset = biomart_output %>% dplyr::rename('ID' = ensembl_gene_id) %>%
               left_join(dataset %>% dplyr::select(ID, contains('MM.')), by='ID')

  
  ##############################################################################
  # PERFORM ENRICHMENT
  
  # Following https://yulab-smu.github.io/clusterProfiler-book/chapter8.html
  
  modules = dataset$Module[dataset$Module!='gray'] %>% as.character %>% table %>% names
  nPerm = 1e5 # 100 times more than the default
  
  enrichment_GO = list()         # Gene Ontology
  enrichment_DO = list()         # Disease Ontology
  enrichment_DGN = list()        # Disease Gene Networks
  enrichment_KEGG = list()       # Kyoto Encyclopedia of Genes and Genomes
  enrichment_Reactome = list()   # Reactome: Pathway db
  
  
  for(module in modules){
    cat('\n')
    cat(paste0('Module: ', which(modules == module), '/', length(modules)))
    geneList = EA_dataset[,paste0('MM.',substring(module,2))]
    names(geneList) = EA_dataset[,'entrezgene'] %>% as.character
    geneList = sort(geneList, decreasing = TRUE)
    
    enrichment_GO[[module]] = gseGO(geneList, OrgDb = org.Hs.eg.db, pAdjustMethod = 'bonferroni', 
                                    pvalueCutoff = 0.1, nPerm = nPerm, verbose = FALSE, seed = TRUE) %>% 
                              data.frame
    enrichment_DO[[module]] = gseDO(geneList, pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1,
                                    nPerm = nPerm, verbose = FALSE, seed = TRUE) %>% data.frame
    enrichment_DGN[[module]] = gseDGN(geneList, pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1,
                                      nPerm = nPerm, verbose = FALSE, seed = TRUE) %>% data.frame
    enrichment_KEGG[[module]] = gseKEGG(geneList, organism = 'human', pAdjustMethod = 'bonferroni', 
                                        pvalueCutoff = 0.1, nPerm = nPerm, verbose = FALSE, seed = TRUE) %>% 
                                data.frame
    enrichment_Reactome[[module]] = gsePathway(geneList, organism = 'human', pAdjustMethod = 'bonferroni', 
                                               pvalueCutoff = 0.1, nPerm = nPerm, verbose = F, seed = T) %>% 
                                    data.frame
    
    # Temporal save, just in case SFARI Genes enrichment fails
    save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, file=file_name)
  }
  

  ##############################################################################
  # PERFROM ENRICHMENT FOR SFARI GENES
  
  # BUILD MAPPING BETWEEN GENES AND SFARI

  # Build TERM2GENE: dataframe of 2 columns with term and gene
  term2gene = biomart_output %>% 
              left_join(genes_info %>% dplyr::select(ID, `gene-score`), 
                         by = c('ensembl_gene_id'='ID')) %>% dplyr::select(-ensembl_gene_id) %>% 
              mutate('SFARI' = ifelse(`gene-score` != 'Others','SFARI','Others'),
                     `gene-score` = ifelse(`gene-score` != 'Others', 
                                           paste0('SFARI Score ',`gene-score`), 'Others')) %>%
              melt(id.vars = 'entrezgene') %>% dplyr::select(value, entrezgene) %>% 
              dplyr::rename('term' = value, 'gene' = entrezgene) %>% distinct
  
  
  # PERFORM GSEA
  enrichment_SFARI = list()
  cat('\n\nGSEA OF SFARI GENES\n')
  
  for(module in modules){
    cat('\n')
    cat(paste0('Module: ', which(modules == module), '/', length(modules)))
    geneList = EA_dataset[,paste0('MM.',substring(module,2))]
    names(geneList) = EA_dataset[,'entrezgene'] %>% as.character
    geneList = sort(geneList, decreasing = TRUE)
      
    enrichment_SFARI[[module]] = clusterProfiler::GSEA(geneList, pAdjustMethod = 'bonferroni',  nPerm = nPerm,
                                                       TERM2GENE = term2gene, pvalueCutoff=1, maxGSSize=2e3,
                                                        verbose = FALSE, seed = TRUE) %>% data.frame
    
    # Temporal save
    save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, 
         enrichment_SFARI, file=file_name)
  }
  
  ##############################################################################
  # PERFROM ENRICHMENT FOR OLD SFARI GENES
  
  # BUILD MAPPING BETWEEN GENES AND SFARI

  # Build TERM2GENE: dataframe of 2 columns with term and gene
  term2gene = biomart_output %>% 
              left_join(genes_info %>% dplyr::select(ID, `old-gene-score`), 
                         by = c('ensembl_gene_id'='ID')) %>% dplyr::select(-ensembl_gene_id) %>% 
              mutate('SFARI' = ifelse(`old-gene-score` != 'Others','SFARI','Others'),
                     `old-gene-score` = ifelse(`old-gene-score` != 'Others', 
                                           paste0('SFARI Score ',`old-gene-score`), 'Others')) %>%
              melt(id.vars = 'entrezgene') %>% dplyr::select(value, entrezgene) %>% 
              dplyr::rename('term' = value, 'gene' = entrezgene) %>% distinct
  
  
  # PERFORM GSEA
  enrichment_old_SFARI = list()
  cat('\n\nGSEA OF OLD SFARI GENES\n')
  
  for(module in modules){
    cat('\n')
    cat(paste0('Module: ', which(modules == module), '/', length(modules)))
    geneList = EA_dataset[,paste0('MM.',substring(module,2))]
    names(geneList) = EA_dataset[,'entrezgene'] %>% as.character
    geneList = sort(geneList, decreasing = TRUE)
      
    enrichment_old_SFARI[[module]] = clusterProfiler::GSEA(geneList, pAdjustMethod = 'bonferroni',  
                                                           nPerm = nPerm, TERM2GENE = term2gene, pvalueCutoff=1, 
                                                           maxGSSize=2e3, verbose = FALSE, seed = TRUE) %>% 
                                     data.frame
    
    # Temporal save
    save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, 
         enrichment_SFARI, enrichment_old_SFARI, file=file_name)
  }
  ##############################################################################
  # Save enrichment results
  save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, 
       enrichment_SFARI, enrichment_old_SFARI, file=file_name)
  
  rm(getinfo, mart, biomart_output, module, term2gene, geneList, EA_dataset, nPerm)
}

# Rename lists
GSEA_GO = enrichment_GO
GSEA_DGN = enrichment_DGN
GSEA_DO = enrichment_DO
GSEA_KEGG = enrichment_KEGG
GSEA_Reactome = enrichment_Reactome
GSEA_SFARI = enrichment_SFARI
GSEA_old_SFARI = enrichment_old_SFARI

file_name = './../Data/ORA.RData'

if(file.exists(file_name)){
  load(file_name)
} else {
  
  ##############################################################################
  # PREPARE DATASET
  
  # Create dataset with top modules membership and removing the genes without an assigned module
  EA_dataset = data.frame('ensembl_gene_id' = genes_info$ID, module = genes_info$Module)  %>% 
               filter(genes_info$Module!='gray')
  
  # Assign Entrez Gene Id to each gene
  getinfo = c('ensembl_gene_id','entrezgene')
  mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl', 
                 host='feb2014.archive.ensembl.org')
  biomart_output = getBM(attributes=getinfo, filters=c('ensembl_gene_id'), 
                         values=genes_info$ID[genes_info$Module!='gray'], mart=mart)
  
  EA_dataset = biomart_output %>% dplyr::rename('ID' = ensembl_gene_id) %>%
               left_join(dataset %>% dplyr::select(ID, Module), by='ID')

  
  ##############################################################################
  # PERFORM ENRICHMENT
  
  # Following https://yulab-smu.github.io/clusterProfiler-book/chapter8.html
  
  modules = dataset$Module[dataset$Module!='gray'] %>% as.character %>% table %>% names
  universe = EA_dataset$entrezgene %>% as.character
  
  enrichment_GO = list()         # Gene Ontology
  enrichment_DO = list()         # Disease Ontology
  enrichment_DGN = list()        # Disease Gene Networks
  enrichment_KEGG = list()       # Kyoto Encyclopedia of Genes and Genomes
  enrichment_Reactome = list()   # Reactome: Pathway db
  
  
  for(module in modules){
    
    genes_in_module = EA_dataset$entrezgene[EA_dataset$Module == module]
    
    enrichment_GO[[module]] = enrichGO(gene = genes_in_module, universe = universe, OrgDb = org.Hs.eg.db, 
                                       ont = 'All', pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1, 
                                       qvalueCutoff = 1) %>% data.frame
    enrichment_DO[[module]] = enrichDO(gene = genes_in_module, universe = universe, qvalueCutoff = 1,
                                       pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1) %>% data.frame
    enrichment_DGN[[module]] = enrichDGN(gene = genes_in_module, universe = universe, qvalueCutoff = 1,
                                         pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1) %>% data.frame
    enrichment_KEGG[[module]] = enrichKEGG(gene = genes_in_module, universe = universe, qvalueCutoff = 1,
                                           pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1) %>% data.frame
    enrichment_Reactome[[module]] = enrichPathway(gene = genes_in_module, universe = universe, qvalueCutoff = 1,
                                                  pAdjustMethod = 'bonferroni', pvalueCutoff = 0.1) %>% 
                                    data.frame
  }
  
  # Temporal save, just in case SFARI Genes enrichment fails
  save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, file=file_name)
  
  
  ##############################################################################
  # PERFROM ENRICHMENT FOR SFARI GENES
  
  # BUILD MAPPING BETWEEN GENES AND SFARI

  # Build TERM2GENE: dataframe of 2 columns with term and gene
  term2gene = biomart_output %>% 
              left_join(genes_info %>% dplyr::select(ID, `gene-score`), 
                         by = c('ensembl_gene_id'='ID')) %>% dplyr::select(-ensembl_gene_id) %>% 
              mutate('SFARI' = ifelse(`gene-score` != 'Others','SFARI','Others'),
                     `gene-score` = ifelse(`gene-score` != 'Others', 
                                           paste0('SFARI Score ',`gene-score`), 'Others')) %>%
              melt(id.vars = 'entrezgene') %>% dplyr::select(value, entrezgene) %>% 
              dplyr::rename('term' = value, 'gene' = entrezgene) %>% distinct
  
  
  # PERFORM GSEA
  enrichment_SFARI = list()
  
  for(module in modules){
      genes_in_module = EA_dataset$entrezgene[EA_dataset$Module == module]
      
      enrichment_SFARI[[module]] = enricher(gene = genes_in_module, universe = universe, 
                                            pAdjustMethod = 'bonferroni', TERM2GENE = term2gene, 
                                            pvalueCutoff = 1, qvalueCutoff = 1, maxGSSize = 50000) %>% 
                                    data.frame %>% dplyr::select(-geneID,-Description)
  }

  ##############################################################################
  # PERFROM ENRICHMENT FOR SFARI GENES
  
  # BUILD MAPPING BETWEEN GENES AND SFARI

  # Build TERM2GENE: dataframe of 2 columns with term and gene
  term2gene = biomart_output %>% 
              left_join(genes_info %>% dplyr::select(ID, `old-gene-score`), 
                         by = c('ensembl_gene_id'='ID')) %>% dplyr::select(-ensembl_gene_id) %>% 
              mutate('SFARI' = ifelse(`old-gene-score` != 'Others','SFARI','Others'),
                     `old-gene-score` = ifelse(`old-gene-score` != 'Others', 
                                           paste0('SFARI Score ',`old-gene-score`), 'Others')) %>%
              melt(id.vars = 'entrezgene') %>% dplyr::select(value, entrezgene) %>% 
              dplyr::rename('term' = value, 'gene' = entrezgene) %>% distinct
  
  
  # PERFORM GSEA
  enrichment_old_SFARI = list()
  
  for(module in modules){
      genes_in_module = EA_dataset$entrezgene[EA_dataset$Module == module]
      
      enrichment_old_SFARI[[module]] = enricher(gene = genes_in_module, universe = universe, 
                                                pAdjustMethod = 'bonferroni', TERM2GENE = term2gene, 
                                                pvalueCutoff = 1, qvalueCutoff = 1, maxGSSize = 5e4) %>% 
                                       data.frame %>% dplyr::select(-geneID,-Description)
  }
  
  ##############################################################################
  # Save enrichment results
  save(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, 
       enrichment_SFARI, enrichment_old_SFARI, file=file_name)
  
  rm(getinfo, mart, biomart_output, gene, module, term2gene, genes_in_module, EA_dataset, universe, file_name)
}

# Rename lists
ORA_GO = enrichment_GO
ORA_DGN = enrichment_DGN
ORA_DO = enrichment_DO
ORA_KEGG = enrichment_KEGG
ORA_Reactome = enrichment_Reactome
ORA_SFARI = enrichment_SFARI
ORA_old_SFARI = enrichment_old_SFARI

rm(enrichment_GO, enrichment_DO, enrichment_DGN, enrichment_KEGG, enrichment_Reactome, enrichment_SFARI, 
   enrichment_old_SFARI)

Both GSEA and ORA are commonly used to study enrichment in sets of genes, but when using them for studying our modules both have shortcomings:

• GSEA takes into consideration some ordering of the genes, in this case given by their Module Membership, which is correlated to the membership of genes to the module, but has two problems:

• Being a continuous scale, it doesn’t separate by a threshold the genes that are truly in the cluster from the rest

• The Module Membership metric is correlated to the real membership of the module, but this correlation is not perfect: a high MM doesn’t always mean the gene belongs to that module, for example, selecting a random module, in the plot below we can see the MM distribution of genes belonging to that module against the rest of the genes belonging to other modules and, although in general, genes belonging to that module have a higher distribution of MM, there is still a big overlap between the two groups, making MM a noisy metric for performing GSEA

modules = dataset$Module[dataset$Module!='gray'] %>% as.character %>% table %>% names

module = sample(modules,1)

plot_data = dataset %>% dplyr::select(Module, paste0('MM.',gsub('#','',module))) %>% 
            mutate(in_module = substring(Module,2) == gsub('#','',module), selected_module = module) %>%
            mutate(alpha = ifelse(in_module, 0.8, 0.1))
colnames(plot_data)[2] = 'MM'

p = plot_data %>% ggplot(aes(selected_module, MM, color = in_module)) + geom_jitter(alpha = plot_data$alpha) + 
    xlab('') + ylab('Module Membership') + coord_flip() + theme_minimal() + 
    theme(legend.position = 'none')

ggExtra::ggMarginal(p, type = 'density', groupColour = TRUE, groupFill = TRUE, margins = 'x', size=1)

rm(modules, module, p, plot_data)

• ORA doesn’t rank the genes within the cluster, but considers them all together, so if you have an enrichment only in the most important genes in the module, it will probably be diluted by the rest of the not so important genes and ORA won’t be able to identify it

So perhaps it could be useful to use both methods together, since they seem to complement each other’s shortcomings very well, performing the enrichment using both methods and identifying the terms that are found to be enriched by both

Note: Since the enrichment in both methods is quite a stric restriction, we decide to relax the corrected p-value threshold (using Bonferroni correction) to 0.1.

compare_methods = function(GSEA_list, ORA_list){
  
  for(top_module in top_modules){
  
    cat(paste0('  \n  \nEnrichments for Module ', top_module, ' (MTcor=',
               round(dataset$MTcor[dataset$Module==top_module][1],2), '):  \n  \n'))
    
    GSEA = GSEA_list[[top_module]]
    ORA = ORA_list[[top_module]]
    
    cat(paste0('GSEA has ', nrow(GSEA), ' enriched terms  \n'))
    cat(paste0('ORA has  ', nrow(ORA), ' enriched terms  \n'))
    cat(paste0(sum(ORA$ID %in% GSEA$ID), ' terms are enriched in both methods  \n  \n'))

    plot_data = GSEA %>% mutate(pval_GSEA = p.adjust) %>% dplyr::select(ID, Description, NES, pval_GSEA) %>%
                inner_join(ORA %>% mutate(pval_ORA = p.adjust) %>% 
                           dplyr::select(ID, pval_ORA, GeneRatio, qvalue), by = 'ID') 
    
    if(nrow(plot_data)>0){
      print(plot_data %>% mutate(pval_mean = pval_ORA + pval_GSEA) %>% 
                          arrange(pval_mean) %>% dplyr::select(-pval_mean) %>% 
            kable %>% kable_styling(full_width = F))
    }
  } 
}


plot_results = function(GSEA_list, ORA_list){
  
  l = htmltools::tagList()

  for(i in 1:length(top_modules)){
    
    GSEA = GSEA_list[[top_modules[i]]]
    ORA = ORA_list[[top_modules[i]]]
    
    plot_data = GSEA %>% mutate(pval_GSEA = p.adjust) %>% dplyr::select(ID, Description, NES, pval_GSEA) %>%
                inner_join(ORA %>% mutate(pval_ORA = p.adjust) %>% dplyr::select(ID, pval_ORA), by = 'ID')
    
    if(nrow(plot_data)>5){
      min_val = min(min(plot_data$pval_GSEA), min(plot_data$pval_ORA))
      max_val = max(max(max(plot_data$pval_GSEA), max(plot_data$pval_ORA)),0.05)
      ggp = ggplotly(plot_data %>% ggplot(aes(pval_GSEA, pval_ORA, color = NES)) + 
                     geom_point(aes(id = Description)) + 
                     geom_vline(xintercept = 0.05, color = 'gray', linetype = 'dotted') + 
                     geom_hline(yintercept = 0.05, color = 'gray', linetype = 'dotted') + 
                     ggtitle(paste0('Enriched terms in common for Module ', top_modules[i])) +
                     scale_x_continuous(limits = c(min_val, max_val)) + 
                     scale_y_continuous(limits = c(min_val, max_val)) + 
                     xlab('Corrected p-value for GSEA') + ylab('Corrected p-value for ORA') +
                     scale_colour_viridis(direction = -1) + theme_minimal() + coord_fixed())
      l[[i]] = ggp
    }
  }
  
  return(l)
}

KEGG

compare_methods(GSEA_KEGG, ORA_KEGG)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 34 enriched terms
ORA has 3 enriched terms
3 terms are enriched in both methods

ID	Description	NES	pval_GSEA	pval_ORA	GeneRatio	qvalue
hsa05322	Systemic lupus erythematosus	2.383671	0.0076646	0.0000729	15/196	3.48e-05
hsa05203	Viral carcinogenesis	2.039923	0.0081770	0.0000695	20/196	3.48e-05
hsa05034	Alcoholism	2.512285	0.0080476	0.0002834	18/196	9.02e-05

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 18 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods

Reactome

compare_methods(GSEA_Reactome, ORA_Reactome)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 94 enriched terms
ORA has 73 enriched terms
64 terms are enriched in both methods

ID	Description	NES	pval_GSEA	pval_ORA	GeneRatio	qvalue
R-HSA-171306	Packaging Of Telomere Ends	2.797288	0.0288531	0.0000000	14/248	0.0000000
R-HSA-110330	Recognition and association of DNA glycosylase with site containing an affected purine	2.809859	0.0289917	0.0000000	14/248	0.0000000
R-HSA-110331	Cleavage of the damaged purine	2.809859	0.0289917	0.0000000	14/248	0.0000000
R-HSA-73927	Depurination	2.809859	0.0289917	0.0000000	14/248	0.0000000
R-HSA-110328	Recognition and association of DNA glycosylase with site containing an affected pyrimidine	2.690112	0.0291747	0.0000000	14/248	0.0000000
R-HSA-110329	Cleavage of the damaged pyrimidine	2.690112	0.0291747	0.0000000	14/248	0.0000000
R-HSA-73928	Depyrimidination	2.690112	0.0291747	0.0000000	14/248	0.0000000
R-HSA-3214842	HDMs demethylate histones	2.459638	0.0292524	0.0000000	16/248	0.0000000
R-HSA-73929	Base-Excision Repair, AP Site Formation	2.684661	0.0292524	0.0000000	14/248	0.0000000
R-HSA-73728	RNA Polymerase I Promoter Opening	2.946047	0.0293359	0.0000000	14/248	0.0000000
R-HSA-5334118	DNA methylation	2.918911	0.0294813	0.0000001	14/248	0.0000000
R-HSA-427359	SIRT1 negatively regulates rRNA expression	2.759346	0.0295267	0.0000001	14/248	0.0000000
R-HSA-5625886	Activated PKN1 stimulates transcription of AR (androgen receptor) regulated genes KLK2 and KLK3	2.863767	0.0295749	0.0000002	14/248	0.0000000
R-HSA-5693571	Nonhomologous End-Joining (NHEJ)	2.303682	0.0295919	0.0000000	15/248	0.0000000
R-HSA-606279	Deposition of new CENPA-containing nucleosomes at the centromere	2.682341	0.0295919	0.0000002	14/248	0.0000000
R-HSA-774815	Nucleosome assembly	2.682341	0.0295919	0.0000002	14/248	0.0000000
R-HSA-1221632	Meiotic synapsis	2.611073	0.0296705	0.0000003	14/248	0.0000000
R-HSA-2299718	Condensation of Prophase Chromosomes	2.945113	0.0297640	0.0000005	14/248	0.0000000
R-HSA-212300	PRC2 methylates histones and DNA	2.784749	0.0297640	0.0000005	14/248	0.0000000
R-HSA-2559586	DNA Damage/Telomere Stress Induced Senescence	2.500765	0.0298261	0.0000001	15/248	0.0000000
R-HSA-5693565	Recruitment and ATM-mediated phosphorylation of repair and signaling proteins at DNA double strand breaks	2.401013	0.0298393	0.0000001	15/248	0.0000000
R-HSA-157579	Telomere Maintenance	2.304934	0.0298393	0.0000014	14/248	0.0000000
R-HSA-5693606	DNA Double Strand Break Response	2.400206	0.0298997	0.0000002	15/248	0.0000000
R-HSA-427389	ERCC6 (CSB) and EHMT2 (G9a) positively regulate rRNA expression	2.815317	0.0298997	0.0000002	15/248	0.0000000
R-HSA-912446	Meiotic recombination	2.545407	0.0300142	0.0000001	16/248	0.0000000
R-HSA-3214858	RMTs methylate histone arginines	2.283929	0.0304344	0.0000000	17/248	0.0000000
R-HSA-73772	RNA Polymerase I Promoter Escape	2.484752	0.0304344	0.0000185	14/248	0.0000004
R-HSA-5250924	B-WICH complex positively regulates rRNA expression	2.530772	0.0304344	0.0000185	14/248	0.0000004
R-HSA-1912408	Pre-NOTCH Transcription and Translation	2.398050	0.0304712	0.0000003	16/248	0.0000000
R-HSA-201722	Formation of the beta-catenin:TCF transactivating complex	2.398704	0.0304712	0.0000026	15/248	0.0000001
R-HSA-8936459	RUNX1 regulates genes involved in megakaryocyte differentiation and platelet function	2.299876	0.0305792	0.0000039	15/248	0.0000001
R-HSA-5625740	RHO GTPases activate PKNs	2.429326	0.0306756	0.0000069	15/248	0.0000002
R-HSA-3214815	HDACs deacetylate histones	2.660197	0.0306822	0.0000006	16/248	0.0000000
R-HSA-977225	Amyloid fiber formation	2.458757	0.0308546	0.0000021	16/248	0.0000001
R-HSA-1912422	Pre-NOTCH Expression and Processing	2.327395	0.0309407	0.0000052	16/248	0.0000001
R-HSA-1500620	Meiosis	2.194344	0.0309407	0.0000052	16/248	0.0000001
R-HSA-5578749	Transcriptional regulation by small RNAs	2.122131	0.0308383	0.0001460	14/248	0.0000025
R-HSA-5250913	Positive epigenetic regulation of rRNA expression	2.450169	0.0309586	0.0000331	15/248	0.0000006
R-HSA-73854	RNA Polymerase I Promoter Clearance	2.387812	0.0309990	0.0000456	15/248	0.0000009
R-HSA-427413	NoRC negatively regulates rRNA expression	2.362381	0.0308939	0.0001704	14/248	0.0000028
R-HSA-2559582	Senescence-Associated Secretory Phenotype (SASP)	2.574914	0.0311734	0.0000011	17/248	0.0000000
R-HSA-73864	RNA Polymerase I Transcription	2.386248	0.0311425	0.0000534	15/248	0.0000010
R-HSA-5250941	Negative epigenetic regulation of rRNA expression	2.314879	0.0309407	0.0002675	14/248	0.0000044
R-HSA-5617472	Activation of anterior HOX genes in hindbrain development during early embryogenesis	2.300994	0.0312880	0.0000167	16/248	0.0000004
R-HSA-5619507	Activation of HOX genes during differentiation	2.300994	0.0312880	0.0000167	16/248	0.0000004
R-HSA-2559580	Oxidative Stress Induced Senescence	2.402468	0.0314089	0.0000308	16/248	0.0000006
R-HSA-1474165	Reproduction	2.096910	0.0314604	0.0000552	16/248	0.0000010
R-HSA-211000	Gene Silencing by RNA	1.888874	0.0315485	0.0003883	15/248	0.0000060
R-HSA-9018519	Estrogen-dependent gene expression	2.322324	0.0319318	0.0000882	17/248	0.0000016
R-HSA-3214847	HATs acetylate histones	2.025256	0.0319366	0.0003048	16/248	0.0000049
R-HSA-68875	Mitotic Prophase	2.185369	0.0319366	0.0016379	15/248	0.0000234
R-HSA-8878171	Transcriptional regulation by RUNX1	1.808700	0.0340641	0.0003342	21/248	0.0000053
R-HSA-212165	Epigenetic regulation of gene expression	2.129608	0.0319318	0.0025336	15/248	0.0000350
R-HSA-2559583	Cellular Senescence	2.065143	0.0331844	0.0014131	18/248	0.0000205
R-HSA-8939211	ESR-mediated signaling	2.398064	0.0337253	0.0017596	19/248	0.0000247
R-HSA-9006931	Signaling by Nuclear Receptors	2.100650	0.0346960	0.0251258	19/248	0.0003138
R-HSA-73884	Base Excision Repair	2.101587	0.0603541	0.0000103	14/248	0.0000002
R-HSA-69473	G2/M DNA damage checkpoint	2.059731	0.0609424	0.0000026	15/248	0.0000001
R-HSA-73886	Chromosome Maintenance	2.028023	0.0616677	0.0001063	14/248	0.0000018
R-HSA-8939236	RUNX1 regulates transcription of genes involved in differentiation of HSCs	1.859668	0.0629809	0.0004414	15/248	0.0000066
R-HSA-157118	Signaling by NOTCH	1.763084	0.0683002	0.0001010	22/248	0.0000018
R-HSA-3247509	Chromatin modifying enzymes	1.678713	0.0696326	0.0133438	20/248	0.0001739
R-HSA-4839726	Chromatin organization	1.678713	0.0696326	0.0133438	20/248	0.0001739
R-HSA-3214841	PKMTs methylate histone lysines	1.986687	0.0897702	0.0000029	14/248	0.0000001

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 27 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods

Gene Ontology

compare_methods(GSEA_GO, ORA_GO)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 31 enriched terms
ORA has 65 enriched terms
0 terms are enriched in both methods

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 29 enriched terms
ORA has 1 enriched terms
0 terms are enriched in both methods

Disease Ontology

compare_methods(GSEA_DO, ORA_DO)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 29 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 18 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods

Disease Gene Network

compare_methods(GSEA_DGN, ORA_DGN)

Enrichments for Module #79AF00 (MTcor=0.87):

GSEA has 6 enriched terms
ORA has 8 enriched terms
4 terms are enriched in both methods

ID	Description	NES	pval_GSEA	pval_ORA	GeneRatio	qvalue
umls:C0278504	Non-small cell lung cancer stage I	2.110432	0.0732876	0.00e+00	16/346	0.0e+00
umls:C0278484	Malignant neoplasm of colon stage IV	2.180269	0.0739020	1.00e-07	16/346	0.0e+00
umls:C0206659	Embryonal Carcinoma	1.955250	0.0789265	3.00e-06	20/346	6.0e-07
umls:C0036421	Systemic Scleroderma	1.668825	0.0973237	5.65e-05	36/346	8.8e-06

Enrichments for Module #00B813 (MTcor=-0.91):

GSEA has 0 enriched terms
ORA has 0 enriched terms
0 terms are enriched in both methods

Session info

sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] WGCNA_1.69             fastcluster_1.1.25     dynamicTreeCut_1.63-1 
##  [4] org.Hs.eg.db_3.8.2     AnnotationDbi_1.46.1   IRanges_2.18.3        
##  [7] S4Vectors_0.22.1       Biobase_2.44.0         BiocGenerics_0.30.0   
## [10] DOSE_3.10.2            ReactomePA_1.28.0      clusterProfiler_3.12.0
## [13] biomaRt_2.40.5         kableExtra_1.1.0       knitr_1.28            
## [16] doParallel_1.0.15      iterators_1.0.12       foreach_1.5.0         
## [19] polycor_0.7-10         expss_0.10.2           ggpubr_0.2.5          
## [22] magrittr_1.5           GGally_1.5.0           gridExtra_2.3         
## [25] viridis_0.5.1          viridisLite_0.3.0      RColorBrewer_1.1-2    
## [28] dendextend_1.13.4      plotly_4.9.2           glue_1.4.1            
## [31] reshape2_1.4.4         forcats_0.5.0          stringr_1.4.0         
## [34] dplyr_1.0.0            purrr_0.3.4            readr_1.3.1           
## [37] tidyr_1.1.0            tibble_3.0.1           ggplot2_3.3.2         
## [40] tidyverse_1.3.0       
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_1.1.0      RSQLite_2.2.0         htmlwidgets_1.5.1    
##   [4] grid_3.6.3            BiocParallel_1.18.1   munsell_0.5.0        
##   [7] codetools_0.2-16      preprocessCore_1.46.0 miniUI_0.1.1.1       
##  [10] withr_2.2.0           colorspace_1.4-1      GOSemSim_2.10.0      
##  [13] highr_0.8             rstudioapi_0.11       ggsignif_0.6.0       
##  [16] labeling_0.3          urltools_1.7.3        polyclip_1.10-0      
##  [19] bit64_0.9-7           farver_2.0.3          vctrs_0.3.1          
##  [22] generics_0.0.2        xfun_0.12             R6_2.4.1             
##  [25] graphlayouts_0.7.0    bitops_1.0-6          reshape_0.8.8        
##  [28] fgsea_1.10.1          gridGraphics_0.5-0    assertthat_0.2.1     
##  [31] promises_1.1.0        scales_1.1.1          ggraph_2.0.3         
##  [34] nnet_7.3-14           enrichplot_1.4.0      ggExtra_0.9          
##  [37] gtable_0.3.0          tidygraph_1.2.0       rlang_0.4.6          
##  [40] splines_3.6.3         lazyeval_0.2.2        acepack_1.4.1        
##  [43] impute_1.58.0         broom_0.5.5           europepmc_0.4        
##  [46] checkmate_2.0.0       BiocManager_1.30.10   yaml_2.2.1           
##  [49] modelr_0.1.6          crosstalk_1.1.0.1     backports_1.1.8      
##  [52] httpuv_1.5.2          qvalue_2.16.0         Hmisc_4.4-0          
##  [55] tools_3.6.3           ggplotify_0.0.5       ellipsis_0.3.1       
##  [58] ggridges_0.5.2        Rcpp_1.0.4.6          plyr_1.8.6           
##  [61] base64enc_0.1-3       progress_1.2.2        RCurl_1.98-1.2       
##  [64] prettyunits_1.1.1     rpart_4.1-15          cowplot_1.0.0        
##  [67] haven_2.2.0           ggrepel_0.8.2         cluster_2.1.0        
##  [70] fs_1.4.0              data.table_1.12.8     DO.db_2.9            
##  [73] triebeard_0.3.0       reprex_0.3.0          reactome.db_1.68.0   
##  [76] matrixStats_0.56.0    xtable_1.8-4          mime_0.9             
##  [79] hms_0.5.3             evaluate_0.14         XML_3.99-0.3         
##  [82] jpeg_0.1-8.1          readxl_1.3.1          compiler_3.6.3       
##  [85] crayon_1.3.4          htmltools_0.4.0       later_1.0.0          
##  [88] Formula_1.2-3         lubridate_1.7.4       DBI_1.1.0            
##  [91] tweenr_1.0.1          dbplyr_1.4.2          MASS_7.3-51.6        
##  [94] rappdirs_0.3.1        Matrix_1.2-18         cli_2.0.2            
##  [97] igraph_1.2.5          pkgconfig_2.0.3       rvcheck_0.1.8        
## [100] foreign_0.8-76        xml2_1.2.5            webshot_0.5.2        
## [103] rvest_0.3.5           digest_0.6.25         graph_1.62.0         
## [106] rmarkdown_2.1         cellranger_1.1.0      fastmatch_1.1-0      
## [109] htmlTable_1.13.3      curl_4.3              shiny_1.4.0.2        
## [112] graphite_1.30.0       lifecycle_0.2.0       nlme_3.1-147         
## [115] jsonlite_1.7.0        fansi_0.4.1           pillar_1.4.4         
## [118] lattice_0.20-41       fastmap_1.0.1         httr_1.4.1           
## [121] survival_3.1-12       GO.db_3.8.2           UpSetR_1.4.0         
## [124] png_0.1-7             bit_1.1-15.2          ggforce_0.3.1        
## [127] stringi_1.4.6         blob_1.2.1            latticeExtra_0.6-29  
## [130] memoise_1.1.0